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1  INTRODUCTION 

This  demonstration  is  designed  to  illustrate  the  discrimination  performance  at  a  challenging  live- 
site  of  a  suite  of  advanced  electromagnetic  induction  (EMI)  modeling  approaches  developed  to 
go  beyond  the  simple  dipole  model  in  accuracy  and  predictive  ability.  The  core  of  the  suite 
consists  of  the  orthonormalized  volume  magnetic  source  (ONVMS)  model  for  target 
characterization,  a  target-counting  pre-processing  procedure  based  on  joint  diagonalization  (JD), 
and  an  implementation  of  the  differential  evolution  (DE)  algorithm  for  nonlinear  optimization 
used  to  locate  targets.  The  study  used  cued  data  sets  collected  at  Camp  Butner  in  North  Carolina 
using  two  next-generation  EMI  sensors,  the  Geometries  MetalMapper  (MM)  and  the  Time- 
domain  Electro-Magnetic  Towed  Array  Detection  System  (TEMTADS)  developed  by  the  NRL 
and  G&G  Sciences.  The  site  was  contaminated  with  fuzes  and  a  mix  of  37-mm  and  105-mm 
munitions.  Each  data  set  was  inverted  with  the  purpose  of  estimating  the  number  of  targets 
producing  each  anomaly  and  the  parameters  associated  with  each  target,  both  extrinsic — its 
orientation,  location  and  depth — and  intrinsic — its  total  volume  magnetic  source  amplitude 
(ONVMS),  which  depends  on  its  size,  shape  and  material  properties.  The  inverted  intrinsic 
parameters  were  then  used  to  classify  the  targets,  and  in  the  end  we  generated  sensor-specific 
dig-lists  for  each  EMI  instrument  and  submitted  them  to  the  Institute  of  Defense  Analyses  (IDA) 
for  independent  scoring. 

1.1  Background 

The  Environmental  Security  Technology  Certification  Program  (ESTCP)  recently  launched  a 
series  of  live-site  UXO  blind  tests  taking  place  in  increasingly  challenging  and  complex  sites  [1- 
4],  The  first  classification  study  was  conducted  in  2007  at  the  UXO  live-site  at  the  former  Camp 
Sibert  in  Alabama  using  two  commercially  available  first- generation  EMI  sensors  (the  EM61- 
MK2  and  the  EM-63,  both  from  Geonics)  [1].  At  that  site,  the  discrimination  test  was  relatively 
simple:  one  had  to  discriminate  large  intact  4.2”  mortars  from  smaller  range  scrap,  shrapnel  and 
cultural  debris,  and  the  anomalies  were  very  well  separated. 

The  second  ESTCP  discrimination  study  took  place  in  2009  at  the  live-UXO  site  at  Camp  San 
Luis  Obispo  (SLO)  in  California  and  featured  a  more  challenging  topography  and  a  wider  mix  of 
targets  of  interest  (TOI)  [4].  Magnetometers  and  first-generation  EMI  sensors  (again  the  Geonics 
EM61-MK2)  were  deployed  on  the  site  and  used  in  survey  mode  for  a  first  screening. 
Afterwards,  two  advanced  EMI  sensing  systems — the  Berkeley  UXO  Discriminator  (BUD)  and 
the  Naval  Research  Laboratory’s  TEMTADS  array — were  used  to  perform  cued  interrogation  of 
a  number  of  the  anomalies  detected.  A  third  advanced  system,  the  Geometries  MetalMapper,  was 
used  in  both  survey  and  cued  modes  for  anomaly  identification  and  classification.  Among  the 
munitions  buried  at  SLO  were  60-mm,  81-mm,  and  4.2"  mortars  and  2.36”  rockets;  three 
additional  types  of  munitions  were  discovered  during  the  course  of  the  demonstration. 

The  third  site  was  chosen  to  be  the  former  Camp  Butner  in  North  Carolina  [2],  This 
demonstration  was  designed  to  investigate  evolving  classification  methodologies  at  a  site  densely 
contaminated  with  small  UXO  (in  this  case  mostly  37-mm  projectiles). 
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1.2  Brief  site  history 

Please  refer  to  the  ESTCP  Live  Site  Demonstration  Plan  [3]. 

1.3  Objective  of  the  demonstration 

The  advanced  EMI  models  we  present  here  (ONVMS  and  JD)  were  developed  under  SERDP 
Project  MM-1572  [5]  and  tested  against  TEMTADS  data  sets  collected  at  the  Aberdeen  Proving 
Ground  (APG)  test  site  in  Maryland  and  in  a  retrospective  analysis  of  cued  TEMTADS  data 
collected  at  the  San  Luis  Obispo  site  [6-10].  The  present  test  of  discrimination  performance 
considers  data  taken  at  Camp  Butner;  this  live  site  was  densely  contaminated  with  small  targets 
such  as  37 -mm  projectiles  and  fuzes,  adding  another  level  of  complexity  into  the  classification 
and  thus  further  demonstrating  the  robustness  of  the  advanced  EMI  models  for  live-site  UXO 
discrimination. 

Overall,  the  principal  objective  of  this  demonstration  was  to  demonstrate  the  models’ 
classification  performance  for  live-site  UXO  problems.  The  specific  technical  objectives  were  to: 

1.  Demonstrate  the  advanced  EMI  models’  classification  accuracy  and  their  applicability  to 
live-site  UXO  discrimination  problems. 

2.  Illustrate  and  document  the  robustness  of  the  data  inversion  and  discrimination  models. 

3.  Invert  targets’  intrinsic  parameters  and  identify  robust  classification  features. 

4.  Indentify  all  seeded  and  native  UXO. 

5.  Eliminate  at  least  75%  of  the  targets  that  do  not  correspond  to  TOI. 

6.  Indentify  sources  of  uncertainty  in  the  classification  process  and  include  them  in  a  dig/no¬ 
dig  decision  process. 

7.  Understand  and  document  the  applicability  and  limitations  of  the  advanced  EMI 
discrimination  technologies  in  the  context  of  project  objectives,  site  characteristics,  and 
suspected  ordnance  contamination. 
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2  TECHNOLOGY 

The  advanced  EMI  models  and  statistical  signal  processing  approaches  developed  and  tested 
over  the  past  three  years  under  SERDP  Project  MM- 1572  [5]  were  able  to  detect  and  identify 
buried  UXO  ranging  in  caliber  from  25  mm  up  to  155  mm.  The  technique  was  seen  to  be 
physically  complete,  fast,  accurate,  and  clutter-tolerant,  and  provided  excellent  classification  in 
both  single-  and  multiple-target  scenarios  when  combined  with  multi-axis/transmitter/receiver 
sensors  like  TEMTADS  and  MetalMapper  [11].  We  start  our  technology  description  with  an 
overview  in  Section  2. 1  of  the  orthonormalized  volume  magnetic  source  model,  which  we  use  to 
represent  the  signatures  and  extract  the  properties  of  multiple  subsurface  targets  simultaneously 
in  an  efficient  manner.  We  then  discuss  in  Section  2.2  a  data  pre-processing  approach  based  on 
joint  diagonalization  that  often  allows  one  to  make  certain  judgments  on  the  number  and  type  of 
unknown  targets  without  the  need  to  perform  inversion.  In  Section  2.3  we  discuss  our  complete 
inversion  procedure,  which  combines  ONVMS  with  Differential  Evolution  optimization,  and  we 
conclude  in  Section  2.4  by  presenting  our  classification  procedures. 

2.1  The  orthonormalized  volume  magnetic  source  model 

The  advanced  models  we  have  developed  for  UXO  discrimination  include  the  normalized 
surface  magnetic  source  (NSMS)  model  [12]  and  the  orthonormalized  volume  magnetic  source 
(ONVMS)  model  [13].  The  ONVMS  model  can  be  considered  as  a  generalized  volume  dipole 
model:  in  it,  an  object’s  response  to  a  sensor  is  modeled  mathematically  using  a  set  of  equivalent 
point-like  analytic  solutions  of  the  Maxwell  equations  (usually  dipoles,  though  charges  are  also  a 
possibility)  distributed  over  a  computational  volume  located  under  an  EMI  sensor  and  potentially 
containing  anomalies.  The  amplitudes  of  the  sources  are  proportional  to  the  component  of  the 
primary  magnetic  field;  once  this  dependence  is  normalized  out,  the  ONVMS  strengths  are 
determined  directly  from  the  data  using  a  set  of  orthogonal  functions. 

Overall,  we  make  the  usual  EMI  assumptions:  we  neglect  displacement  currents  everywhere,  as 
well  as  electric  fields  and  conduction  currents  in  air  and  soil.  The  primary  magnetic  field 
generated  by  the  sensor  penetrates  the  objects  in  its  vicinity  to  some  degree,  inducing  eddy 
currents  and  magnetic  dipoles  inside  them,  which,  in  turn,  produce  a  secondary  or  scattered 
magnetic  field.  This  is  the  field  that  we  propose  to  represent  as  being  due  to  a  volumetric 
distribution  of  magnetic  dipoles: 

Hsc  (r,  P)  =  j  — !-3  (3RR  -  T)  m(r; ,  p)dv'  =  J  G(r,  <)  •  m(rv',  p)dV ,  ( 1 ) 

V  y 

where  p  =  {t,f}  is  time  or  frequency,  R  is  the  unit  vector  along  R  =  r-rj,  r(  is  the  position 
of  the  v'  -th  infinitesimal  dipole  in  the  volume  V,  r  is  the  observation  point,  and  I  and  G(r,r') 
are  respectively  the  identity  and  Green  dyads.  The  induced  magnetic  dipole  moment  m(rj,/r)  at 
point  rj  on  the  surface  is  related  to  the  primary  field  through  m(r',  p)  =  M(r',p)  •  Hpr(r') ,  where 
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M(rJ \,p)  is  the  symmetric  polarizability  tensor.  The  secondary  magnetic  field  at  any  point  can 
be  expanded  into  a  set  of  orthonormal  functions  ^.(r)  as 

Nr 

H(r)  =  X^(R)-b.,  (2) 

(=i 

where  we  have  also  introduced  the  expansion  coefficients  b. .  The  \j/  are  linear  combinations  of 

dipole  Green  dyads  and  are  guaranteed  to  be  orthonormal  as  a  result  of  the  Gram-Schmidt 
orthogonalization  procedure;  thanks  to  this  property,  the  amplitudes  of  the  tensor  elements 

M.(p)  can  be  determined  without  having  to  solve  a  linear  system  of  equations.  Two  great 

advantages  of  ONVMS  are  that  it  takes  into  account  the  mutual  couplings  between  different 
sections  of  the  targets  and  that  it  avoids  matrix  singularity  problems  in  multi-object  cases.  It 
treats  single-  and  multi-target  scenarios  on  the  same  footing.  Once  the  tensor  elements  and 
locations  of  the  responding  dipoles  are  determined,  one  can  group  them  in  space  and  determine 
the  total  polarizability  tensor  within  each  group,  which  is  then  jointly-diagonalized  in  time  to 
extract  the  temporal  decay  law  of  its  diagonal  elements.  These  diagonal  elements  have  been 
shown  to  be  intrinsic  to  the  objects,  and  can  be  used,  either  on  their  own  or  in  combination  with 
other  quantities,  in  discrimination  processing  [14].  The  theoretical  basis  of  the  ONVMS  model  is 
outlined  in  Appendix  C. 

2.2  Joint  diagonalization  data  preprocessing 

Advanced  electromagnetic  induction  (EMI)  sensors  currently  feature  multi-axis  illumination  of 
targets  and  tri-axial  vector  sensing,  or  exploit  multi-static  array  data  acquisition  [11].  They 
produce  data  of  high  density,  quality,  and  diversity,  and  have  been  combined  with  advanced  EMI 
models  to  provide  superb  classification  performance  [14]  relative  to  the  previous  generation  of 
single-axis  monostatic  sensors  [15-17].  To  take  advantage  of  the  rich  data  sets  that  these  sensors 
provide,  we  recently  developed  and  successfully  demonstrated  a  discrimination-oriented  data 
pre-processing  scheme  based  on  joint  diagonalization  (JD)  [13].  Let  us  illustrate  the  method  by 
describing  its  TEMTADS  implementation.  TEMTADS  consists  of  25  transmit/receive  pairs  of 
square  coil  antennas,  each  consisting  of  a  35-centimeter  (cm)  transmitter  loop  surrounding  a 
concentric  25-cm  receiver  coil,  arranged  in  a  5x5  square  grid.  The  sensor  activates  the 
transmitter  loops  in  sequence,  one  at  a  time,  and  for  each  transmitter  all  receivers  receive, 
measuring  the  complete  transient  response  over  a  wide  dynamic  range  of  time  ranging  from 
approximately  100  microseconds  (ps)  to  25  milliseconds  (ms)  and  distributed  over  Nq  =121  time 
channels.  The  sensor  thus  provides  25  x  25  spatial  data  points  at  any  given  time  channel  tq, 
q  =  1,2,...,  Nq.  If  we  define  Hk,m  as  the  z-component  of  the  magnetic  field  measured  by  the  m-th 
receiver  coil  when  the  k- th  transmitter  is  active,  then  each  row  of  the  measured  multi-static 
response  (MSR)  data  matrix 
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is  a  measured  data  vector  for  the  k- th  transmitter,  where  k  =  1,2,  . . K.  For  the  TEMTADS 
system  M  =  25  is  the  number  of  receivers  and  K  =  25  is  the  number  of  transmitters.  For  each  time 
channel  the  M  x  K  MSR  matrix  can  be  decomposed  into  its  eigenvectors  U(f9)  and  eigenvalues 
D(t<?)  using  the  singular  value  decomposition  (SVD)  to  obtain 


H(g  =  U(t?)D(gUr(g.  (4) 

However,  in  order  to  relate  the  eigenvalues  to  the  number  of  potential  targets,  we  need  to  find  a 

N 

unitary  matrix  V  of  eigenvectors  shared  by  all  [H(t  )}  ’  matrices  that  simultaneously  removes 
all  their  off-diagonal  elements: 


D(t9)  =  VrH(gV,  q  =  1,  . . Nq.  (5) 

In  general,  it  is  not  the  case  that  the  matrix  V  will  cancel  all  the  off-diagonal  elements  of  all  the 
D(tq),  but  a  unitary  V  can  be  sought  that  minimizes  the  sum  of  their  squares;  this  is  the  gist  of  the 
JD  approach  [18].  The  diagonal  elements  of  D(f9)  are  the  time-dependent  eigenvalues  of  the 
measured  MSR  matrix  and  contain  information  about  the  targets  that  contribute  to  a  given  signal. 
Our  studies  show  that  three  diagonal  elements  of  the  MSR  matrix  usually  suffice  to  describe  one 
target.  We  have  also  found  that  JD  is  a  robust  technique  for  extracting  signals  due  to  targets  for 
data  with  low  signal  to  noise  ratios.  See  Appendix  C  for  a  more  detailed  exposition  of  the 
method. 

2.3  EMI  Data  inversion:  A  global  optimization  technique 

Determining  a  buried  object’s  orientation  and  location  is  a  non-linear  problem.  Inverse-scattering 
problems  are  solved  by  defining  an  objective  function  [5]  that  measures  the  mismatch  between 
modeled  and  measured  magnetic-field  data  and  proceeding  to  find  its  global  minimum.  Standard 
gradient  search  approaches  often  suffer  from  a  profusion  of  local  minima  that  sometimes  result 
in  incorrect  location  and  orientation  estimates.  To  avoid  this  problem  we  recently  adapted  and 
incorporated  a  different  class  of  global  optimization  search  algorithms,  among  them  differential 
evolution  (DE)  [19-20],  a  heuristic,  parallel,  direct-search  method  for  minimizing  non-linear 
functions  of  continuous  variables  that  is  straightforward  to  implement  and  has  good  convergence 
properties.  We  have  combined  DE  with  ON  VMS  to  invert  digital  geophysical  EMI  data  [5,  14]. 
All  EMI  optimizations  are  split  into  linear  and  nonlinear  parts,  alternating  between  the  two  and 
iterating  to  minimize  the  objective  function.  Once  the  target  locations  are  found,  the  amplitudes 
of  the  responding  ONVMS  are  determined  and  used  to  classify  the  object  relative  to  the  targets 
of  interest. 
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2.3.1  Discrimination  parameters 

To  classify  targets  in  this  demonstration  we  used  ONVMS  combined  with  DE  optimization  and 
JD  to  invert  for  the  locations  and  electromagnetic  signatures  of  the  TOI.  The  model  provides  at 
least  three  independent  principal-axis  total  ONVMS  parameters  per  target  that  can  be  used  for 
discrimination.  The  total  time-dependent  ONVMS  depends  on  the  size,  geometry,  and  material 
composition  of  the  object  in  question:  1)  Early  time  gates  bring  out  the  high-frequency  response 
to  the  shutdown  of  the  exciting  field;  the  induced  eddy  currents  in  this  range  are  superficial,  and 
a  large  total  ONVMS  amplitude  at  early  times  correlates  with  large  objects  and  large  surface 
area.  2)  At  late  times,  when  the  eddy  currents  have  diffused  completely  into  the  object  and  low- 
frequency  harmonics  dominate,  the  EMI  response  relates  to  the  metal  content  (i.e.,  the  volume) 
of  the  target.  Thus  a  smaller  but  compact  target  has  a  relatively  weak  early  response  that  dies 
down  slowly,  while  a  large  but  thin  or  hollow  object  has  a  strong  initial  response  that  decays 
quickly.  Thus  the  extracted  ONVMS  parameters  can  be  used  to  form  feature  vectors  for 
classification. 

2.3.2  Classification  approaches 

The  power-law/exponential-decay  parameters  extracted  from  total  ONVMS  time-decay  curves 
tend  to  follow  definite  patterns  when  TOI  of  the  same  kind  are  interrogated  under  different 
conditions.  The  parameters  thus  tend  to  cluster  together  in  ways  that  provide  clues  as  to  the 
features  of  the  different  TOI  present  in  a  survey,  and  by  comparing  total  ONVMS  parameters  of 
unknown  objects  to  those  of  previously  characterized  targets  one  can  predict  the  class/cluster  to 
which  the  unknown  targets  belong. 

There  are  many  clustering  techniques  available,  such  as  K-means  [21],  principal  component 
analysis  [22],  and  support  vector  machines  [18,  21,  23-28],  which  are  largely  heuristically 
motivated  and  do  not  require  an  underlying  statistical  model.  A  possible  alternative  is  the  model- 
based  clustering  method,  which  is  based  on  the  assumption  that  the  data  are  composed  of  a  finite 
mixture  of  different  distributions  of  the  same  type  (e.g.,  multivariate  Gaussians)  but 
characterized  by  different  sets  of  parameters.  This  procedure  has  the  obvious  advantage  that  it  is 
possible  to  choose  the  optimal  number  of  clusters  and  the  distributions  that  best  fit  the  data  using 
some  objective  statistical  criterion — such  as  the  Akaike  Information  Criterion  (AIC)  or  the 
Bayesian  Information  Criterion  (BIC) — while  for  the  other  methods  those  questions  are  open  to 
discussion. 

Clustering  methods  can  also  be  categorized  into  unsupervised  and  supervised.  Supervised 
clustering  uses  parameters  extracted  from  a  set  of  training/calibration  samples,  whereas 
unsupervised  clustering  applies  the  same  classification  criteria  to  all  targets,  regardless  of  size, 
composition,  and  decay  curves.  The  former  is  obviously  advantageous  in  that  it  utilizes  as  prior 
knowledge  additional  information  from  the  training  data.  Several  supervised  clustering 
techniques  like  support  vector  machines  [18,  21,  23-28]  and  template  matching  have  been 
applied  to  UXO  discrimination.  A  straightforward  approach  to  implement  model-based 
supervised  clustering  is  to  (a)  estimate  the  parameters  from  the  training  sample  and  (b)  use  the 
estimated  values  of  the  parameters  to  classify  the  TOI  in  the  blind  test  dataset.  A  drawback  of 
this  method  is  that  only  the  training  set  is  used  to  estimate  the  parameters  and  information  from 
the  blind  dataset  is  completely  ignored.  To  avoid  this,  during  this  study  a  semi-supervised 
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classification  algorithm  was  applied  in  order  to  create  custom  training  sets.  Decay  parameters 
such  as  Maa  (A),  and  Maa (t\ )  /  Maa(tn),  where  t„  is  a  given  time  channel  and  Maa  is  the  total 
ONVMS  along  the  «-th  axis,  were  used  to  cluster  the  inverted  parameters  into  classes.  We 
assumed  that  there  were  K  clusters  and  that  each  of  them  was  described  by  a  parametric 
continuous  or  discrete  distribution  (e.g.,  a  Gaussian).  This  let  us  arrange  total  ONVMS-extracted 
time-decay  parameters  in  the  n  x  m  matrix  Y  =  [Yi,  Y2,  . . .,  Y m],  with  n  the  number  of  anomalies 
and  m  the  number  of  chosen  parameters.  Taking  each  Y,-  to  follow  an  m-dimensional  mixture  of 
normal  distributions  allowed  us  to  express  the  total  ^-cluster  distribution  as 

F<Y,)  =  ZV,<Y,  (6) 

k=\ 

where  wk,  is  the  mixing  weight  of  cluster  k,  S;A=|  vv,  =  1 ,  and 


(7) 


is  the  probability  density  of  the  k-th  normal  distribution  with  a  mean  vector  \ik  (an  m  x  1  vector) 
and  a  variance-covariance  matrix  ak  (an  m  x  m  matrix).  The  mixing  weight  wk  is  defined  as  the 
proportion  of  anomalies  that  belong  to  the  k-th  cluster.  The  parameters  \ik,  <sk,  and  wk  are 
estimated  by  the  maximum  likelihood  (ML)  criterion  using  the  expectation  maximization 
algorithm  [29]. 


2.3.3  Classification  using  template  matching 


Template  matching  is  a  classification  approach  that  identifies  an  unknown  target  by  comparing 
its  extracted  features — here  the  total  ONYMS — to  those  of  a  set  stored  in  a  reference  library.  The 
comparison  can  be  carried  out  either  1)  by  using  code  that  computes  least-squares  mismatches  or 
2)  by  visual  inspection.  To  avoid  false  negatives  and  classify  targets  accurately  we  used  both 
approaches. 


2.4  Details  of  classification  schemes 

The  discrimination  process  comprises  three  sequential  tasks:  data  collection,  data  inversion,  and 
classification.  Each  EMI  sensor  produces  unique  data  sets  and  therefore  requires  its  own  data 
inversion  and  classification  schemes.  This  section  summarizes  the  data  inversion  and 
classification  schemes  for  the  5  x  5  TEMTADS  array  and  the  MM  sensor. 
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Figure  1.  TEMTADS  multi-static  response  matrix  eigenvalues  versus  time  for  (top  row)  a 
105-mm  HE  projectile  and  a  105-mm  HEAT  round,  (center  row)  an  M-48  fuze  and  a  37-mm 
munition,  and  (third  row)  two  clutter  scenarios,  one  with  two  items  (left)  and  another  with 
several  (right). 
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Figure  2.  TEMTADS  multi-static  response  matrix  eigenvalues  versus  time  for  some  samples  of 
requested  anomalies. 


2.4.1  Camp  Butner  TEMTADS  data  inversion  and  classification  scheme 

The  modeling  approach  used  to  invert  TEMTADS  data  with  the  ONVMS-DE  algorithm  is 

described  in  detail  in  [5].  Here  we  summarize  the  main  steps  employed  to  invert  and  classify 

TEMTADS  data  for  the  particular  case  of  Camp  Butner. 

Step  1.  Data  pre-processing:  All  TEMTADS-data  were  pre-processed  using  a  Matlab  code  (see 
Appendix  C  in  [30]  that  reads  comma-delimited  CSV  files  and  translates  them  to  ASCII  files 
compatible  with  the  ONVMS-DE  code  (ONVMSJVIM.exe).  The  user  needs  only  specify  the 
path  to  the  folder  with  the  CSV  files;  the  code  then  converts  them  all. 

Step  2.  Use  equation  (3)  to  construct  the  TEMTADS  MSR  matrix  H(r9). 

Step  3.  Eigenvalue  analysis:  The  JD  algorithm  constructs  a  multi-static  response  matrix  using 
TEMTADS  data  and  computes  its  eigenvectors  and  eigenvalues,  the  latter  as  a  function  of 
time;  these  are  depicted  for  some  of  the  Camp  Butner  anomalies  in  Figure  1  and  Figure  2. 
The  eigenvalues  of  the  MSR  data  matrix  are  intrinsic  properties  of  the  targets,  and  each  target 
has  at  least  three  eigenvalues  above  the  threshold  (the  noise  level  is  composed  of  low- 
magnitude  eigenvalues).  For  example,  Figure  1  shows  the  eigenvalues  extracted  for  a  105- 
mm  HE  projectile,  a  105-mm  HEAT  round,  an  M-48  fuze,  a  37-mm  projectile,  and  some 
clutter  items.  Each  target  is  seen  to  have  distinguishable  eigenvalues,  and  we  used  these  to 
make  an  initial  classification.  Note  that  the  magnitudes  of  the  MSR  eigenvalues  depend  on 
the  depths  and  orientations  of  the  targets  [5],  so  the  user  must  rely  only  on  their  shapes  when 
performing  JD-based  classification.  As  the  number  of  targets  increases  (as  in  Figure  2  and 
the  third  row  of  Figure  1),  so  does  the  number  of  above-noise  eigenvalues.  We  examined  the 
eigenvalues’  time-decay  curves  for  each  case  and  used  them  to  estimate  the  corresponding 
number  of  targets. 

Step  4.  Once  we  estimated  the  number  of  targets  and  SNR  for  each  anomaly  we  inverted  all  cued 
MM  datasets  using  a  multi-target  combined  ONVMS-DE  algorithm.  This  gave  us  the 
extrinsic  and  intrinsic  parameters  for  all  targets,  including  the  total  ONVMS  as  shown  in 
Figure  3  and  Figure  4. 
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Step  5.  Create  a  custom  training  list  using  eigenvalue  time-decay  curves  and  request  the  ground 
truth.  For  the  most  part,  the  JD-based  list  contained  either  the  anomalies  that  had  too  many 
above-threshold  eigenvalues,  like  the  samples  depicted  in  the  third  row  of  Figure  1  and  in 
Figure  2,  or  anomalies  which  had  very  small  eigenvalues.  We  requested  two  batches  of 
training  data.  The  first  batch  contained  65  anomalies,  all  of  which  were  clutter;  some  had  six 
eigenvalues  above  the  noise  level,  while  others  had  several  eigenvalues  mixed  with  the  noise. 
The  second  batch  consisted  of  10  anomalies,  of  which  8  corresponded  to  UXO. 
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Figure  3:  Inverted  total  ONVMS  time-decay  profiles  for  four  Camp  Butner  targets:  (top  row) 
105-mm  HE  munition  and  105-mm  HEAT  round,  and  (bottom)  M-48  Fuze  and  37-mm  projectile 
with  copper  band. 
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Figure  4:  Inverted  total  ONVMS  time-decay  profiles  for  37-mm  projectiles  without  copper  band. 


Step  6.  Once  we  had  the  ground  truth  for  all  75  custom  identified  anomalies  we  proceeded  to 
classify  all  TEMTADS  targets  using  the  inverted  total  ONVMS  as  discriminating  features. 

Step  7.  Create  ranked  dig  list.  Armed  with  the  custom  identified  training  list  and  the  inverted 
total  ONVMS  for  each  case  we  created  a  library  for  M-48  fuzes  and  37-mm  projectiles 
without  copper  band.  We  did  not  request  training  data  for  either  of  the  105-mm  UXO  or  for 
the  37-mm  projectile  with  copper  band  because  we  already  had  TEMTADS  test-stand  data 
for  these  targets.  The  inverted  total  ONVMS  for  the  anomalies  that  were  classified  as  TOI 
appear  in  Figure  3  and  Figure  4.  All  the  inverted  total  ONVMS  are  seen  to  cluster  well,  and 
each  target  has  a  total  ONVMS  with  features — such  as  its  amplitude  at  the  first  time  channel, 
its  decay  rate,  or  the  separation  between  the  (blue)  primary  and  (red  and  green)  secondary 
components  at  different  time  channels — that  make  it  amenable  to  identification.  (The  most 
difficult  differences  to  discern  were  between  the  M-48  fuzes  of  Figure  3  and  the  37-mm 
projectiles  without  copper  band  of  Figure  4).  These  features  allowed  us  to  classify  targets  as 
UXO  or  clutter  and  also  let  us  sort  the  UXO  by  caliber.  With  this  knowledge  we  created  a 
prioritized  dig  list  that  we  cross-validated  using  the  time-decay  curves  of  the  JD  eigenvalues. 

Step  8.  Submit  the  dig  list  to  ESTCP.  The  final  prioritized  dig  list  was  submitted  to  the  Institute 
for  Defense  Analyses  (IDA)  for  independent  scoring.  The  scored  results  were  sent  back  in 
the  form  of  a  receiver  operating  characteristic  (ROC)  curve,  which  appears  in  Figure  5.  We 
can  see  that  a)  of  the  75  targets  that  were  dug  for  training,  68  targets  were  not  TOI  (shift 
along  x-axis)  and  seven  were  UXO  (shift  along  y-axis);  b)  for  95%  TOI  classification  (the 
pink  dot  in  Figure  5)  only  seven  extra  (false  positive)  digs  are  needed;  c)  to  classify  all  TOI 
correctly  (the  light  blue  dot)  only  21  extra  (false  positive)  digs  are  needed;  d)  in  order  to 
increase  the  classification  confidence,  the  algorithm  requested  an  additional  thirty  digs  after 
all  TOI  had  been  identified  correctly. 
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Figure  5.  ROC  curve  for  the  Camp  Butner  TEMTADS  blind  test. 


Difficult  TOIs: 
#460  (37mm) 
#2504  (37mm) 
#2405  (37mm) 
#2340  (37mm) 
#2023  (37mm) 
#2017  (37mm) 
#1981  (37mm) 
#1945  (37mm) 
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Figure  6:  Left :  Scatter  plot  for  all  MM  anomalies  based  on  the  extracted  total  ONVMS.  Right: 
Probability  function  for  all  MM  anomalies. 

2.4.2  Camp  Butner  MM  data  inversion  and  classification  scheme 

The  Geometries  MetalMapper  (MM)  is  a  commercially  available  advanced  EMI  sensor.  Its 
geometric  configuration  is  described  in  [11],  and  the  ONVMS-DE-based  modeling  approach  for 
inverting  the  data  it  provides  is  described  in  detail  in  [5].  The  Camp  Butner  MM  cued  data  sets 
were  collected  by  teams  from  Sky  Research  and  Geometries.  A  data  analysis  revealed  that  the 
Geometries  MM  data  had  a  smaller  SNR  than  the  Sky  Research  data.  In  addition,  it  was  found 
that  the  y-receiver  of  sensor  #3  was  not  working  properly  when  Geometries  took  their 
measurements  [11].  To  accommodate  this  prior  information  it  was  necessary  to  exclude  this 
receiver  from  all  Geometries  MM  data  sets  and  perform  two  kinds  of  classification  (statistical 
and  library-matching)  to  provide  cross-checks.  Moreover,  it  was  necessary  to  pay  extra  attention 
to  the  Geometries  MM  data  when  performing  library  (fingerprint)  matching:  any  anomaly  whose 
primary  total  ONVMS  resembled  the  total  primary  ONVMS  of  any  TOI  was  ranked  as  a  TOI  or 
included  in  the  custom  training  list.  Below  we  summarize  the  main  steps  taken  during  inversion 
and  classification  of  Camp  Butner  MM  data. 

Step  1.  Extract  the  total  ONVMS  for  each  anomaly:  We  ran  the  program  ONVMS_MM.EXE  to 
extract  the  parameters  for  all  targets. 

Step  2.  Create  a  custom  training  list:  For  classification  we  used  the  ratio  of  the  inverted  total 
ONVMS  at  the  30th  time  channel  to  that  at  the  first.  The  values  of  logio  [M::(t\ )  /  AT-tCo)] 
versus  logio  [Mzz(ti)]  are  plotted  on  Figure  6  (left)  for  all  the  data.  They  clearly  exhibit  a  wide 
spread  of  values.  To  use  these  features  for  statistical  classification,  and  for  determining 
clusters  and  a  classification  probability  function,  we  first  divided  the  scatter  plot  of  Figure  6 
(left)  into  1 1  subsections  (shown  in  the  figure).  Then  to  each  of  these  subsections  we  applied 
a  Gaussian  mixture  model  (implemented  as  part  of  Matlab’s  Statistics  Toolbox)  assuming 
that  it  contained  five  clusters.  (The  number  of  subsections  and  the  number  of  clusters  are  not 
critical  at  this  early  stage.  We  selected  the  quoted  values  to  ensure  that  all  clusters  were 
found  and  at  the  same  time  to  minimize  the  number  of  anomalies  whose  ground  truth  was 
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requested.)  From  the  Gaussian  mixture  model  we  found  the  mean  and  standard  deviation  for 
each  cluster  and  built  a  global  classification  probability  function,  depicted  in  Figure  6  (right). 

Step  3.  Request  the  ground  truth  for  selected  anomalies:  Using  the  classification  probability 
function  we  created  a  first  custom  training  dig  list  that  contained  55  anomalies  (i.e.,  one 
anomaly  for  each  cluster  identified  in  the  preceding  step)  and  requested  the  ground  truth.  The 
success  of  classification  depends  on  the  selection  of  features,  the  separation  between 
different  classes  in  feature  space,  and  the  ability  of  the  sensor  data  to  constrain  the  estimated 
features.  We  would  like  to  emphasize  that  in  some  cases,  due  to  poor  signal-to-noise  ratio  (a 
particularly  common  occurrence  for  Camp  Butner  Geometries  MM  data),  the  feature  vectors 
from  UXO  targets  were  corrupted  (see  Figure  7)  or  were  similar  to  those  from  clutter.  In  such 
cases,  and  given  our  previous  studies  (SLO),  we  recognized  that  automatic  statistical 
discrimination  algorithms  had  limitations.  As  a  consequence  we  decided  to  have  the  final 
classification  decision  be  made  by  expert  judgment,  overriding  automated  classification  if 
necessary.  This  was  achieved  by  inverting  all  the  MM  data  using  the  combined  ONVMS-DE 
algorithm  as  though  there  were  one,  two  or  three  targets  present  and  comparing  the  resulting 
total  ONVMS  amplitudes  case-by-case.  Whenever  we  spotted  significant  differences  we 
examined  the  curves  visually,  like  in  the  sample  case  of  Figure  7,  and,  based  on  this 
examination,  requested  the  ground  truth  for  an  additional  60  Geometries  MM  anomalies. 

Create  a  ranked  dig  list:  Using  the  requested  ground  truth  for  121  training  anomalies  and  the 
the  inverted  total  ONVMS  for  each  case  we  created  libraries  for  the  105 -mm  and  37-mm 
projectiles  and  the  M48  fuzes.  The  inverted  total  ONVMS  for  the  anomalies  that  were  correctly 
correctly  classified  as  TOI  appear  in  Figure  8  and  Figure  9.  The  ground  truth  from  the  custom 
training  data  set  also  let  us  classify  all  targets  as  either  TOI  or  not  TOI  using  the  probability 
function  of  Figure  6  (right).  The  classification  based  on  the  supervised  clustering  is  plotted  in 

Step  4.  Figure  10:  the  red  circles  correspond  to  TOI  and  the  green  dots  to  clutter. 

Step  5.  Submit  the  dig  list  to  ESTCP:  Using  the  clustering  and  library-matching  techniques  we 
classified  the  anomalies  as  TOI  or  not  TOI.  The  ranked  list  was  submitted  to  the  IDA  for 
scoring.  The  results  are  shown  on  Figure  11. 
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Figure  7:  Inverted  magnetic  dipole  polarizability  (left)  and  total  ONVMS  (right)  time-decay 
profiles  for  MM  anomaly  #2504.  The  thin  red  lines  show  a  library  sample,  while  the  thick  blue 
and  green  lines  show  the  inversion  results. 
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Figure  8.  Inverted  total  ONVMS  time-decay  profiles  from  Camp  Butner  MM  data  sets:  (top  row) 
105-mm  HE  shells  and  105  mm  HEAT  rounds,  (bottom  row)  M-48  Fuzes  and  37-mm  projectiles. 
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Figure  9.  Inverted  total  ONVMS  time-decay  profiles  from  Camp  Butner  MM  data  sets  for 
37-mm  projectiles  without  driving  copper  bands. 
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Figure  10:  Result  of  the  supervised  clustering  classification  for  the  Camp  Butner  MM  anomalies 
using  the  logarithms  of  Mzz(t\)  /  Mzz(ti{))  and  Mzz(t\).  The  supervised  clustering  was  trained  with 
calibration  data.  The  green  markers  correspond  to  clutter  and  the  red  ones  to  TOI. 
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Butner  Dartmouth  Advanced  Models  None  MetalMapper  Custom  v2  TOI 


Difficult  TOIs: 
#884  (37mm) 
#178  (37mm) 
#1728 (Fuze) 
#1002  (37mm) 
#1000 (Fuze) 
#856  (37mm) 
#773  (37mm) 
#720  (37mm) 


Figure  1 1 :  ROC  curve  for  the  Camp  Butner  MetalMapper  blind  test. 


Step  6.  The  scored  results  for  the  2291  Camp  Butner  MM  anomalies,  depicted  in  Figure  11,  show 
that  a)  of  the  121  targets  that  were  dug  for  training,  120  targets  were  not  TOI  (shift  along  x- 
axis)  and  one  was  (shift  along  y-axis);  b)  for  95%  TOI  classification  (pink  dot  in  Figure  11) 
eight  extra  (false  positive)  digs  are  needed;  c)  to  classify  all  TOI  correctly  (light  blue  dot) 
only  32  additional  digs  are  needed;  d)  for  increased  classification  confidence  the  algorithm 
requested  33  additional  digs  after  all  the  TOI  were  identified  correctly. 
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3  PERFORMANCE  OBJECTIVES 

The  performance  objectives  of  this  ESTCP  live  site  discrimination  study  were:  to  achieve  high 
probability  of  discrimination  of  UXO  from  among  a  wide  spread  of  clutter;  to  process  all  data 
sets;  to  minimize  the  number  of  data  that  could  not  be  analyzed  or  decided  upon;  to  minimize  the 
number  of  false  positives;  and  to  identify  all  UXO  with  high  confidence.  The  performance 
objectives  are  summarized  in  Table  1. 


Table  1:  Performance  objectives 


Performance 

Objective 

Metric 

Data  Required 

Success  Criteria 

Maximize  correct 
classification  of 
munitions 

Number  of  targets  of 
interest  retained 

•  Prioritized  anomaly 
lists 

•  Scoring  reports  from 
the  Institute  for 
Defense  Analyses 
(IDA) 

The  approach  correctly 
classifies  all  targets  of 
interest 

Maximize  correct 
classification  of  non¬ 
munitions 

Number  of  false  alarms 
eliminated 

•  Prioritized  anomaly 
lists 

•  Scoring  reports  from 
the  IDA 

Reduction  of  false  alarms 
by  over  75%  while 
retaining  all  targets  of 
interest 

Specification  of  no-dig 
threshold 

Probability  of  correct 
classification  and 
number  of  false  alarms 
at  demonstrator 
operating  point 

•  Demonstrator- 
specified  threshold 

•  Scoring  reports  from 
the  IDA 

Threshold  specified  by  the 
demonstrator  to  achieve 
the  criteria  specified 
above 

Minimize  the  number 
of  anomalies  that 
cannot  be  analyzed 

Number  of  anomalies 
that  must  be  classified 
as  “Unable  to  Analyze” 

•  Demonstrator  target 
parameters 

Reliable  target  parameters 
can  be  estimated  for  over 
90%  of  anomalies  on  each 
sensor’s  detection  list. 

Correct  estimation  of 
target  parameters 

Accuracy  of  estimated 
target  parameters 

•  Demonstrator  target 
parameters 

•  Results  of  intrusive 
investigation 

Total  ONVMS  ±  10% 

X,  Y  <  +  10  cm 

Z  <  +  5  cm 

size  ±  10% 
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3.1  Objective:  maximize  correct  classification  of  munitions 

An  effective  technology  for  discrimination  of  munitions  should  maximize  the  number  of  targets 
of  interest  (TOI)  it  can  classify  as  such  (thus  distinguishing  them  from  non-TOI)  with  high 
confidence. 

3.1.1  Metric 

Identify  all  seeded  and  native  TOI  with  high  confidence  using  advanced  EMI  discrimination 
technologies.  (The  Program  Office  did  not  quantify  “high  confidence.”  Our  estimates  were  based 
on  using  the  extracted  total  ONVMS  as  input  to  statistical  classification  algorithms  and  expert 
judgment.  Every  anomaly  that  was  close  to  a  TOI  cluster  in  feature  space  and  had  /  >  10“8 , 
where  /  is  the  probability  density  function  of  equation  (7),  was  considered  a  possible  TOI;  the 
expert  then  inspected  the  corresponding  TONVMS  curve  for  symmetry  (manifested  by  equal 
secondary  and  tertiary  TONVMS)  and  signal-to-noise  ratio.) 

3.1.2  Data  requirements 

We  analyzed  data  from  two  instruments,  the  5  x  5  TEMTADS  array  and  the  MetalMapper.  For 
each  sensor  we  identified  custom  training  data  sets.  We  requested  the  ground  truth  for  the  custom 
training  data  sets  and  used  it  to  validate  the  models  for  each  specific  sensor.  We  generated  dig- 
lists  that  were  scored  by  the  IDA. 

3.1.3  Success  criteria  evaluation  and  results 

The  objective  was  considered  to  be  met  if  all  seeded  and  native  UXO  items  could  be  identified 
below  an  analyst-specified  no-dig  threshold. 

3.1.4  Results 

The  objective  was  successfully  met.  All  TOI,  both  seeded  and  native,  were  identified  using  our 
advanced  EMI  discrimination  technology.  The  ROC  curves  for  the  Camp  Butner  test 
demonstrate  that  all  TOI  were  classified  correctly,  both  for  TEMTADS  (Figure  5)  and  the 
MetalMapper  (Figure  11). 

3.2  Objective:  maximize  correct  classification  of  non-munitions 

The  technology  aims  to  minimize  the  number  of  false  negatives,  i.e.,  maximize  the  correct 
classification  of  non-TOI. 

3.2.1  Metric 

We  compared  the  number  of  non-TOI  targets  that  can  be  left  in  ground  with  high  confidence 
using  the  advanced  EMI  discrimination  technology  to  the  total  number  of  false  targets  that  would 
be  present  if  the  technology  were  absent. 
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3.2.2  Data  requirements 

This  objective  required  prioritized  anomaly  lists,  which  our  team  generated  independently  for 
each  sensor,  and  for  its  evaluation  we  needed  scoring  reports  from  the  IDA. 

3.2.3  Success  criteria  evaluation  and  results 

The  objective  was  considered  to  have  been  met  if  the  method  eliminated  at  least  75%  of  targets 
that  did  not  correspond  to  TOI  in  the  discrimination  step. 

3.2.4  Results 

This  objective  was  successfully  met.  The  advanced  forward  EMI  models  were  able  to  extract 
robust  classification  parameters.  Using  the  extracted  parameters  the  classification  algorithms 
were  able  to  eliminate  at  least  93%  of  non-TOI  targets,  as  stated  by  the  IDA. 

3.3  Objective:  specify  a  no-dig  threshold 

This  project  aims  to  provide  a  high-confidence  classification  approach  for  UXO-site  managers. 
A  critical  quantity  for  minimizing  the  residual  risk  of  UXO  and  providing  regulators  with 
acceptable  confidence  is  a  specific  no-dig  threshold. 

3.3.1  Metric 

We  compared  an  analyst’s  no-dig  threshold  point  to  the  point  where  100%  of  munitions  were 
correctly  identified. 

3.3.2  Data  requirements 

To  meet  this  requirement  we  needed  scoring  reports  from  the  IDA. 

3.3.3  Success  criteria  evaluation  and  results 

The  objective  would  be  met  if  a  sensor-specific  dig  list  placed  all  the  TOI  before  the  no-dig  point 
and  if  additional  digs  (false  positives)  were  requested  after  all  TOI  were  identified  correctly. 

3.3.4  Results 

This  objective  was  successfully  met  for  all  data  sets  (see  Figure  5  and 

Figure  10).  The  stop-dig  threshold  was  specified  based  on  the  mismatch  between  the  primary 
total  ONYMS  for  the  TOI  library  items  and  for  the  test  anomalies.  Both  the  magnitudes  and 
shapes  of  the  time-decay  curves  were  used  for  the  final  classification  and  stop-digging  decisions. 
All  uncertain  or  difficult  targets  were  included  in  the  custom  training  lists. 
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3.4  Objective:  minimize  the  number  of  anomalies  that  cannot  be  analyzed 

Some  anomalies  may  not  be  classified,  either  because  of  the  data  are  not  sufficiently 
informative — the  sensor  physically  cannot  provide  the  data  to  support  classification  for  a  given 
target  at  a  given  depth — or  because  the  data  processing  was  inadequate.  The  former  is  a  measure 
of  instrument  performance  for  all  anomalies  for  which  all  data  analysts  converge.  The  latter  is  a 
measure  of  the  quality  of  a  data  analysis  when  a  target  diagnostic  differs  from  those  made  by 
other  analysts. 

3.4.1  Metric 

The  metric  for  this  objective  is  the  number  of  anomalies  that  cannot  be  analyzed  by  a  particular 
method,  and  the  intersection  of  all  anomaly  lists  among  all  analysts. 

3.4.2  Data  requirements 

Each  analyst  submitted  their  anomaly  list.  IDA  scored  all  lists  and  returned  a  list  of  anomalies 
that  could  not  be  analyzed  by  any  analyst  (“cannot  analyze”  or  “failed  classification”). 

3.4.3  Success  criteria  evaluation  and  results 

The  objective  was  met  if  at  least  95%  of  a  set  of  selected  anomalies  could  be  analyzed. 

3.4.4  Results 

This  objective  was  successfully  met.  All  four  data  sets  for  all  anomalies  were  analyzed.  Not  a 
single  anomaly  was  ranked  as  “cannot  analyze.” 

3.5  Objective:  correct  estimation  of  target  parameters 

The  combined  ONVMS-DE  algorithm  provides  intrinsic  and  extrinsic  parameters  for  the 
different  targets.  The  intrinsic  parameters  were  used  for  classification,  while  the  extrinsic 
parameters  (i.e.,  the  target  locations)  were  utilized  for  residual  risk  assessment. 

3.5.1  Metric 

The  classification  results  entirely  depend  on  how  accurately  these  parameters  are  estimated. 

3.5.2  Data  requirements 

To  achieve  this  objective  we  inverted  and  tabulated  the  intrinsic  and  extrinsic  parameters  for  all 
targets.  To  validate  extracted  extrinsic  parameters  we  needed  results  of  intrusive  investigations. 

3.5.3  Success  criteria  evaluation  and  results 

The  objective  was  met  if  the  targets’  intrinsic  parameters  varied  within  +10%,  the  extracted  x-y 
location  within  +10  cm,  and  the  depth  within  +5  cm. 
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3.5.4  Results 

After  learning  the  ground  truth  for  the  extrinsic  parameters  we  evaluated  the  estimates  given  by 
our  procedure.  We  look  at  target  depth  in  more  detail,  as  it  is  the  easiest  parameter  to  compare. 
Figure  12  and  Figure  13  show  (for  the  MetalMapper  and  TEMTADS,  respectively)  the 
distribution  of  depth  errors  (defined  here  by  I  zmodel  -  Zdata  I )  The  MetalMapper  discrepancies 
have  a  mean  of  4.53  cm  and  a  standard  deviation  of  4.53  cm;  for  TEMTADS  the  mean  is 
4.37  cm  and  the  standard  deviation  is  3.76  cm.  The  agreement  between  inverted  and  actual 
values  is  good. 

The  error  in  horizontal  location,  defined  by  ((xmodel -X^)2 +(Tmodel -Tdata)2)1/2 ,  obeys  a 
similar  distribution.  For  the  MetalMapper  we  find  a  mean  of  8.5  cm  and  a  standard  deviation  of 
10.5  cm,  while  for  TEMTADS  the  mean  is  9.3  cm  and  the  standard  deviation  is  1 1.4  cm. 

We  can  see  that  the  discrepancies  in  depth  are  on  average  smaller  than  5  cm,  while  the  average 
horizontal  error  is  smaller  than  10  cm.  Thus  the  objective  was  successfully  met. 


|  Error|  [cm] 

Figure  12.  Histogram  of  depth  errors  (defined  as  |  zmodel  -  Zdata  I)  for  the  set  of  Camp  Butner 
MetalMapper  anomalies.  The  distribution  shown  has  a  mean  of  4.53  cm  and  a  standard  deviation 
of  4.53  cm.  There  is  good  agreement  between  the  estimates  and  the  ground  truth. 
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|Error|  [cm] 


Figure  13.  Histogram  of  depth  errors  for  the  Camp  Butner  TEMTADS  anomalies.  The  mean 
error  here  is  4.37  cm  and  the  standard  deviation  is  3.76  cm.  Again  we  see  acceptable  agreement 
with  the  ground  truth. 


23 


April  2012 


Demonstration  report 


Advanced  EMI  models  for  Camp  Butner 


4  TEST  DESIGN 

The  only  required  test  at  the  Camp  Butner  site  entailed  collecting  target  characterization  training 
data:  Using  a  calibration  pit,  the  data-collection  team  made  a  series  of  static  measurements  of 
example  targets  at  several  depths  and  attitudes  in  order  to  cross-check  models,  confirm  Tx  and 
Rx  polarity  for  the  sensors,  and  characterize  the  so-called  “library  targets.” 

4.1  Site  preparation 

N/A 

4.2  Demonstration  schedule 


Preparation 

Calibration 

Blind  data  set 

Post-survey 

analysis 

Tasks  and  demonstration  stages 

Aug2010 

Sep 

’10 

Oct 

’10 

Nov 

’10 

Dec 

’10 

Jan 

’ll 

Feb 

’ll 

1 .  Invert  all  calibration  data 

X 

2.  Invert  5x5  TEMTADS  data 

X 

3.  Invert  MM  data 

X 

4.  Build  custom  training  data  sets  and  request 
ground  truth  for  TEMTADS 

X 

5.  Build  custom  training  data  sets  and  request 
ground  truth  for  MM 

X 

6.  Redefine  the  MM  classifier  and  request  more 
training  data  if  necessary 

X 

7.  Redefine  the  TEMTADS  target  classifier  and 
request  additional  training  data  if  necessary 

X 

8.  Generate  MM  dig  list  and  submit  to  IDA 

X 

9.  Generate  TEMTADS  dig  list  and  submit  to 
IDA 

X 

10.  Conduct  retrospective  analysis  if  needed 

X 

X 

REPORTING: 

1 1 .  Draft  demonstration  report 

X 

12.  Final  demonstration  report 

X 

Figure  14.  Gantt  chart  showing  a  detailed  schedule  of  the  activities  conducted  at  Camp  Butner. 
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5  DATA  ANALYSIS  PLAN 

We  analyzed  all  cued  data  for  the  MetalMapper  and  TEMTADS  sensors  and  produced 
prioritized  dig  lists  for  independent  scoring. 

5.1  Extracting  target  locations 

Target  locations  were  determined  relative  to  the  sensor  coordinate  system  using  the  differential 
evolution  algorithm.  Object  responses  were  modeled  with  ONVMS.  This  combined  ONVMS-DE 
algorithm  was  run  for  single-  and  multi-target  cases  and  provided  target  locations. 

5.2  Extracting  target  intrinsic  parameters 

5.2.1  Single  targets 

The  combined  ONVMS-DE  algorithm  yields  the  targets’  intrinsic  total  ONVMS,  which  we  used 
for  classification.  The  total  ONVMS  contains  three  moments,  Mxx(t),  Myy(t),  and  Mzz(t),  along  the 
primary  axes  in  the  target’s  own  reference  frame.  These  moments  are  similar  to  simple  dipole 
moment  components  but  carry  more  information,  accounting  for  the  targets’  inherent 
heterogeneities.  The  ONVMS-DE  algorithm  outputs  the  time-decay  curves  of  the  target’s  total 
ONVMS  tensor  M,/4).  The  next  step  is  to  determine  the  time  decay  of  the  primary  components 
of  the  total  ONVMS  in  the  target’s  reference  frame.  While  this  can  be  done  by  standard 
diagonalization  (i.e.,  finding  M (4)  =  V(4)D(4)VT(4),  where  V(4)  contains  the  eigenvectors  of 
M (4),  it  is  more  convenient  to  perform  a  joint  diagonalization,  M (4)  =  VD(4)VT,  where  now  the 
eigenvectors  are  shared  by  all  time  channels;  this  allows  us  to  extract  more  reliable  total 
ONVMS  values  and  reduce  uncertainty.  The  resulting  temporal  decays  of  the  total  principal 
ONVMS  for  the  Camp  Butner  anomalies  are  illustrated  in  Figure  3  and  Figure  4  for  TEMATDS 
and  in  Figure  8  and  Figure  9  for  the  MetalMapper.  The  results  show  that  the  inverted  parameters 
are  clustered  very  well  for  large  targets.  As  the  target  size  decreases  the  magnitude  of  the 
extracted  total  ONVMS  has  a  progressively  larger  spread  (especially  for  MM  data),  but  the 
shapes  of  the  time-decay  curves  remain  the  same,  indicating  that  both  the  shape  and  the 
magnitude  of  the  ONVMS  should  be  considered  trying  to  achieve  accurate  classification. 

5.2.2  Multi-target  cases 

A  similar  approach  is  carried  out  if  more  than  one  subsurface  target  is  expected.  The  DE 
algorithm  now  searches  for  the  locations  and  the  total  ONVMS  of  several  objects.  Such  multi¬ 
target  inversion  is  crucial  in  the  field  for  cases  in  which  a  signal  from  a  UXO  is  mixed  with  EMI 
signals  from  nearby  clutter.  Our  two-target  inversion  code  yields  three  sets  of  location  and  total 
ONVMS  estimates:  one  for  Target  1,  one  for  Target  2,  and  a  combined  estimate  with  Targets  1 
and  2  represented  by  a  single  object.  (In  the  case  of  3 -target  inversion,  seven  sets  of  data  are 
expected:  only  Target  1,  only  Target  2,  only  Target  3,  Targets  1  and  2  as  a  single  object,  Targets 
2  and  3  as  a  single  object,  Targets  1  and  3  as  a  single  object,  and  all  three  targets  acting  as  a 
single  object.  In  the  general  case  of  n  targets  one  expects  n(n  -  1)  +  1  sets  of  ONVMS  curves). 
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5.3  Selection  of  intrinsic  parameters  for  classification 

Most  UXO  are  bodies  of  revolution,  and  consequently  the  two  secondary  polarizability  elements 
are  degenerate.  However,  live-site  UXO  discrimination  studies  have  repeatedly  shown  that  this 
symmetry  can  be  compromised  due  to  low  SNR,  especially  for  small  or  deep  targets.  A  good 
classification  of  object  features  can  then  be  obtained  by  using  only  the  principal  component  of 
the  total  ON  VMS,  Mzz(t).  Furthermore,  to  limit  the  number  of  relevant  features  for  use  in 
classification  we  extract  parameters  exclusively  from  the  main  polarizability  Mzz(t),  both  to 
represent  size  (via  Mzz(t\))  and  wall  thickness  (via  Mzz(ti )  /  Mzz(tn)).  The  interested  reader  is 
referred  to  Section  2.4. 

5.4  Training 

Our  classification  approach  is  based  on  custom  training  data.  At  the  first  stage  of  the  process  we 
used  a  semi-supervised  clustering  technique  to  indentify  potential  site-specific  TOI.  Below  are 
the  basic  steps  performed  during  training  data  selection;  for  more  details  regarding  each  specific 
sensor  see  Section  2.4. 

(a)  The  targets’  intrinsic  features,  Mzz(t\)  and  M--(t\ )  /  Mzz{t„),  were  selected  from  the 
extracted  total  ONVMS.  The  time  channel  n  was  chosen  based  on  feature  separation. 
EMI  data  sets  were  produced  for  all  anomalies  using  both  single-  and  multi-object 
inversions. 

(b)  Initial  clustering  was  performed.  The  ground  truth  was  requested  for  all  targets  whose 
features  were  located  closest  to  the  corresponding  cluster  centroid  and  had  TOI-like 
ONVMS  features. 

(c)  Clusters  containing  at  least  one  TOI  were  identified,  and  a  smaller  domain  was  selected 
within  the  feature  space  for  further  interrogation. 

(d)  Additional  clustering  was  performed  within  the  selected  domain,  and  those  targets  with 
features  closest  to  the  corresponding  cluster  centroids  were  probed  for  ground  truth.  The 
clusters  with  at  least  one  identified  UXO  were  marked  as  suspicious.  The  total  ONVMS 
curves  were  inspected  within  the  selected  domain. 

(e)  All  targets  whose  features  (based  on  multi-object  inversion  and  library  matching)  fell 
inside  any  of  the  suspicious  clusters  were  used  to  train  the  statistical  classifier  and  the 
library-matching  procedure. 

5.5  Classification 

(f)  Probability  density  functions  were  extracted  for  single-  and  multi-target  scenarios. 

(g)  All  of  the  unknown  targets  were  scored  based  on  the  probability  density  functions. 

(h)  Dig  lists  were  produced  for  both  single-  and  multi-object  cases  and  compared  to  each 
other  to  find  similarities  and  differences. 

(i)  All  items  were  further  analyzed  using  library  matching,  and  all  total  ONVMS  time-decay 
curves  were  inspected  visually. 
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(j)  A  set  of  anomalies  were  identified  and  additional  training  data  sets  were  requested.  The 
new  information  was  incorporated  into  the  Gaussian  mixture  model  and  all  items  were  re¬ 
scored. 

(k)  Based  on  the  previous  steps  a  classification  threshold  was  selected  and  a  final  dig  list  was 
produced. 

5.6  Decision  memo 

The  algorithms  used  to  select  training  data  and  to  perform  inversion  and  classification  for  the 
Camp  Butner  test  are  described  in  Section  2.4.  Using  the  inversion,  clustering,  classification  and 
data-requesting  procedures  outlined  above  we  produced  a  ranked  anomaly  list  formatted  as 
specified  by  the  IDA[31]. 
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6  COST  ASSESSMENT 

Time  and  resources  were  tracked  for  each  task  to  assess  the  cost  of  deploying  the  technology  at 
future  live  sites.  A  summary  of  the  time  needed  to  classify  each  anomaly  successfully  appears  in 
Table  2.  This  assessment  does  not  include  actual  data-taking,  as  our  team  did  not  perform  that 
work.  Neither  does  it  include  computer  runtime,  which  usually  takes  place  “in  the  background” 
(i.e.,  concurrently  with  other  tasks,  overnight,  etc.),  nor  the  one-time  cost  involved  in  developing 
the  code.  The  quoted  values  are,  of  course,  averages:  some  anomalies  can  be  identified  almost 
instantaneously,  while  others  require  much  time  and  effort.  The  estimates  stem  from  personal 
observation — the  requirement  to  systematically  keep  track  of  the  time  was  made  official  after  the 
results  had  been  submitted.  Finally,  the  values  given  correspond  to  the  times  needed  by  a 
seasoned  expert.  Analysts  with  less  experience  initially  took  on  average  twice  as  long  to  identify 
and  classify  each  anomaly,  but  quickly  became  much  faster. 


Table  2:  Cost  model  for  advanced  EMI  model  demonstration  at  the  former  Camp  Butner 


Cost  Category 

Description 

Cost 

Pre-processing 

Time  required  to  perform  eigenvalue  extraction, 
check  data  quality,  and  estimate  the  number  of 
potential  anomalies 

15  sec/anomaly 

Parameter  extraction 

Time  required  extract  target  feature  parameters 

15  sec/anomaly 

Classifier  training 

Time  required  to  optimize  classifier  design  and  train 

30  sec/anomaly 

Classification  and 
construction  of  a 
ranked  anomaly  list 

Time  required  to  classify  anomalies  in  the  test  set 
and  construct  the  ranked  anomaly  list 

30  sec/anomaly 

Total 

90  sec/anomaly 
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7  MANAGEMENT  AND  STAFFING 

Figure  15  is  the  organization  chart  for  the  personnel  involved  in  the  demonstration.  Their 
responsibilities  are  as  follows: 

1.  Fridon  Shubitidze  -  Principal  Investigator.  Developed  and  implemented  most  of  the 
preprocessing  and  inversion  routines  used.  Classified  MetalMapper  data  using  a  Gaussian 
mixture  model. 

2.  Irma  Shamatava  -  Sky  Research  Geophysicist.  Participated  in  the  inversion  and 
classification  of  TEMTADS  data. 

3.  Alex  Bijamov  -  Engineer  at  Dartmouth  College.  Classified  TEMTADS  data  via  semi- 
unsupervised  parameter  clustering. 


Fridon  Shubitidze,  PI: 

MM  data  inversion  and  classification 


Mrs.  Irma  Shamatava: 

TEMTADS  inversion  and  classification 

Dr.  Alex  Bijamov,  Dartmouth  College: 
TEMTADS  data  feature-parameter 
clustering  and  classification 

Figure  15:  Project  management  hierarchy. 
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9  APPENDICES 

9.1  Appendix  A:  Health  and  Safety  Plan  (HASP) 

As  this  effort  does  not  involve  field  data  collection,  no  HASP  is  required. 

9.2  Appendix  B:  Points  of  Contact 

Points  of  contact  (POCs)  involved  in  the  demonstration  and  their  contact  information  are 
presented  in  Table  3. 

Table  3:  Points  of  Contact  for  the  advanced  EMI  models’  demonstration 


POINT  OF 
CONTACT 
Name 

ORGANIZATION 

Name 

Address 

Phone 

Fax 

E-mail 

Role  in 
Project 

Dr.  Fridon 
Shubitidze 

Sky  Research,  Inc. 

Tel:  603  643  2876 

Fax:  603-643-5161 
fridon.  shubitidze  @  skvresearch.com 

pi 

Erik  Russell 

Sky  Research,  Inc. 

3  Schoolhouse  Ln, 
Etna,  NH  03750, 
USA 

Tel:  541-552-5197 

Fax:  603-643-5161 
Erik.Russell@skvresearch.com 

Project 

Coordinator 

Dr.  Herb 
Nelson 

ESTCP  Program 
Office, 

ESTCP  Office 

901  North  Stuart  St, 
Suite  303 
Arlington,  VA 
22203-1821 

Tel:  703-696-8726 

Herbert.Nelson@osd.mil 

ESTCP 

Munitions 

Management 

Program 

Manager 
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9.3  Appendix  C: 

9.3.1  The  orthonormalized  volume  magnetic  source  model 

Most  EMI  sensors  are  composed  of  separate  transmitting  and  receiving  coils.  When  the  operator 
activates  the  sensor,  a  current  runs  through  the  transmitter,  resulting  in  the  establishment  of  a 
(“primary”  or  “principal”)  magnetic  field  in  the  surrounding  space  (Figure  16).  According  to  the 
elementary  atomic  model  of  matter,  all  materials  are  composed  of  atoms,  each  with  a  nucleus 
and  a  cloud  of  orbiting  electrons.  The  electrons  cause  circulating  currents  and  form  microscopic 
magnetic  dipoles.  For  most  materials,  in  the  absence  of  an  external  magnetic  field  the  magnetic 
dipoles  of  atoms  have  random  orientations,  resulting  in  no  magnetic  moment.  According  to 
Faraday’s  law,  an  external  time-varying  magnetic  field  induces  eddy  currents  in  conducting 
bodies  by  an  alignment  of  the  magnetic  moments  of  the  “spinning”  electrons  and  a  magnetic 
moment  due  to  a  change  in  the  orbital  motion  of  electrons.  These  currents  and  magnetization  in 
turn  generate  a  (“secondary”  or  “scattered”)  magnetic  field  that  also  varies  with  time  and  induces 
measurable  currents  in  the  receivers.  The  induced  magnetic  dipoles/eddy  currents  are  distributed 
inside  the  object  and  produce  a  magnetic  field  intensity  H  outside.  The  magnetic  field  due  to  the 
i  -th  source  can  then  be  expressed  at  any  observation  point  r  as  the  matrix-vector 
productEquation  Section  (Next) 


H.(r)  =  G  (r)m  , 


(A.l) 


where  m  is  a  1  x  6  dimensional  vector  whose  components  (Mxxj,  Mxyx,  MXZti,  Mxxx,  Myz>i,  M-.j) 
are  the  elements  of  the  target’s  magnetic  polarizability  tensor  M  s  and  [T]  is  the  3  x  6  matrix 
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whose  elements  are  as  follows: 
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When  there  are  several  such  sources,  the  total  field  can  be  expressed  as  a  superposition: 


M 

H(r)  =  £G(r)in.= 

t=l 


G,  G 


m 


(A.3) 


Before  going  further  we  note  that  our  method  takes  as  input  the  (in  principle  unknown)  number 
M  of  radiating  sources.  For  advanced  EMI  sensors  such  as  the  MetalMapper  and  2x2  and  5x5 
TEMTADS  arrays  we  have  developed  a  procedure  based  on  joint  diagonalization,  sketched  in 
Section  9.3.4,  that  estimates  M  starting  from  raw  data  and  with  no  need  for  inversion.  For  other 
sensors  one  may  proceed  by  letting  M  vary  as  part  of  an  optimization  routine. 

The  superposition  (A.3)  can  be  used  to  carry  out  one-  and  multi-object  inversions  starting  from 
data  taken  at  an  ensemble  of  points.  All  the  measured  H  -values — which  can  pertain  to  multiple 
transmitters,  multiple  receivers,  changing  sensor  locations,  and  different  vector  components — are 
strung  together  in  a  one-dimensional  array,  while  the  corresponding  Green  functions  are  stacked 
as  matrix  rows.  The  resulting  composite  G  matrix  can  then  be  (pseudo)inverted  to  find  the 
strengths  of  the  sources.  This  procedure,  which  is  nothing  other  than  the  dipole  model  if  each 
body  is  taken  to  be  represented  by  one  source  only,  works  well  for  one  or  two  sources,  but  for 
larger  numbers  becomes  very  time-consuming  (since  the  Green  matrix  becomes  very  large)  and 
increasingly  ill-posed,  usually  requiring  regularization.  The  ONVMS  method  is  designed  to 
circumvent  these  difficulties. 
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Figure  16:  A  metallic  object  under  the  transmitter.  The  target’s  EMI  response  at  the  receiver  coil 
can  be  calculated  from  the  equivalent  surface  or  volume  magnetic  dipole  moment  dm. 

9.3.2  Orthonormal  Green  functions 

The  method  starts  from  the  realization  that  the  matrix-vector  product  (A.l)  is  valid  at  any 
observation  point  r  and,  in  particular,  at  every  point  r  on  the  plane  surface  delimited  by 
TEMTADS.  If  we  introduce  the  inner  product 

(A,b)=  f  ATBds=  f  ATBds+\  ATBds  +  ...,  (A.4) 

'  '  J  s  Jrx0  Jrx,  ' 

where  the  integral  is  computed  over  the  “sensitive”  surfaces  of  the  sensor  (which  are  contiguous 
in  the  case  of  TEMTADS,  but  not  necessarily  for  other  instruments),  and  if  furthermore  we  can 
find  a  basis  of  Green  functions  orthogonal  under  this  measure, 

M 

H(r,)  =  X'I'A)b,  suchthat  j^k/  =  Fjdjk ,  (A.5) 

7  =  1 

where  S.k  is  a  Kronecker  delta,  then  it  is  possible  to  find  the  source  amplitudes  b  without  costly 
and  ill-conditioned  inversions  simply  by  exploiting  the  sifting  property  of  the  orthogonal  basis: 

M  M 

=  =  E  W/  =  (A.6) 

7=1  7=1 


and  thus 


"('I', .«)■  (A-7) 

which  clearly  does  not  involve  solving  a  linear  system  of  equations;  it  is  necessary  to  invert  only 
the  6  x  6  matrix  Fk .  Moreover,  this  definition  of  the  coefficients  b  guarantees  that  they  are 
“optimal”  in  the  sense  that  the  expansion  (A.5)  yields  the  least  mean-square  error 
(H-2M,'F.b.,H-2M,vP.b.\  [32], 

\  7=1  7  7  7=1  7  7  /  L  J 

To  construct  the  set  of  orthonormal  Green  functions  we  resort  to  a  generalization  of  the  Gram- 
Schmidt  procedure  [33],  Assuming  that  the  Green  matrices  are  linearly  independent — i.e.,  that 
we  cannot  have  a  collection  of  distinctly  located  dipole  sources  combining  to  produce  no 
measurable  field  unless  their  amplitudes  all  vanish — we  define 
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*1  =  ^, 

m- 1 

¥  =G  (A.8) 

m  m  '  2  k  mk 

k= 1 

M- 1 

=gm  -  y 

M  M  Z-J  /:  Mk  ’ 

fc=l 

where  the  6  x  6  matrices  A  obey  A  =  0  for  /  <  k  .  Enforcing  the  orthogonality  relation  (A.5)  is 

equivalent  to  setting  Z'Pn,Gm)  =  FAm  forn  <  m,  and  using  this  relation  twice  in  definition  (A.8) 
we  find  that 


(  ”-1 

A.„ = F;'  i  c ,  -  2Xf*a*  i  .  (A.?) 

V  ^=i  / 

where  the  overlap  integral  C  =(G  ,G  )  . 

At  the  end  of  the  process  it  is  necessary  to  recover  an  expansion  expressed,  like  (A.l),  in  terms 
of  the  actual  Green  functions,  in  part  because  the  functions  V1V  are  orthogonal  (and  defined)  only 

at  points  on  the  receivers,  and  in  part  because  of  the  non-uniqueness  of  the  coefficients  b  due  to 

the  arbitrary  order  in  which  the  G.  enter  the  recursion  (A.8).  To  that  end,  we  express 


=  <A1°) 

k= 1 

and  to  find  the  coefficients  B  we  compare  expansion  (A.  10)  term  by  term  to  the  definition 
(A.8)  and  use  the  rule  that  A.k  =  0  for  j  <k  to  find 


B  =  I,  the  identity, 

mm  J 

B  =  -A  , 

m— 1 

B  =  -V  B.  A  .  for  1  <  q  <  m  -  2, 

mq  Iq  ml  1 

l=q 

in  terms  of  which  we  recover  the  physical  polarizability  elements: 


(A.ll) 


(A.  12) 
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9.3.3  ONVMS  procedure 

With  all  the  pieces  in  place,  we  can  sketch  an  algorithm  to  invert  EMI  data  using  the  ONVMS 
model: 

1)  Given  a  number  of  sources  and  their  tentative  locations,  find  the  Green  tensors  G.  using 
equation  (A.2)  and  compute  the  overlap  integrals  Gmn  using  the  inner  product  (A.4). 

2)  Determine  the  first  normalization  factor,  Fl  =  \GrGt ) ,  and  use  it  to  find  all  the  Gram- 
Schmidt  coefficients  A  with  n  =  1 :  A  =  F~lC '  . 

mn  ml  1  1  m 

3)  Set  m  =  2 ;  compute,  in  sequence, 

a)  The  coefficients  Amn  with  n  =  2,  ,m-l  using  equation  (A.9); 

b)  The  function  'Pm  using  the  expansion  (A.8); 

c)  The  normalization  factor  F  =(x¥  ,¥  / ; 

increase  m  by  1  and  iterate  until  all  sources  have  been  included. 

4)  Once  all  the  A  ,  F  ,  and  T*  are  known,  find  B  using  (A.  11). 

y  mn  m  m  mq  G  x  7 

5)  Use  the  orthonormality  of  the  new  Green  functions  to  determine  the  source  amplitudes 
using  =  F“1^vF^,Hdata\,  as  in  (A.7).  Take  the  measured  field  to  be  piecewise 
constant — i.e.,  constant  throughout  each  receiver — when  evaluating  the  integrals. 

6)  Use  the  computed  b  ,  B  ,  and  Gm,  along  with  the  expansion  (A.  12),  to  generate  the 
secondary  field  prescribed  by  the  given  number  of  sources  at  the  given  locations. 

7)  Compare  the  model  prediction  with  the  measured  data,  vary  the  source  locations,  and 
iterate  until  the  least-squares  discrepancy  between  prediction  and  measurement  attains  a 
suitable  minimum. 

The  procedure  as  written  applies  to  only  one  time  gate,  but  the  extension  to  fully  time-dependent 
functions  is  straightforward:  we  need  only  substitute  the  vectors  b  and  Hdata  for  two- 

dimensional  arrays  where  the  columns  denote  time.  The  relations  between  the  two,  namely  (A.7) 
and  (A.  12),  acquire  multiple  right-hand-sides,  and  the  optimization  mentioned  on  Step  7  of  the 
algorithm  is  constrained  further.  As  a  final  remark  we  note  that  rigorously  speaking  the 
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coefficients  b  (and,  for  that  matter,  the  amplitudes  mt )  are  not  the  polarizabilities  themselves 
but  relate  more  closely  to  their  time  derivatives  [34-35]. 

The  great  advantage  of  the  ONVMS  technique  is  that  it  takes  into  account  mutual  couplings 
between  different  parts  of  targets  and  avoids  matrix  singularity  problems  in  cases  with  multiple 
objects.  Once  the  polarizability  tensor  elements  and  the  locations  of  the  elemental  responding 
dipoles  are  determined  one  can  group  them  according  to  their  volume  distribution.  For  each 
group  a  total  polarizability  tensor  can  be  computed  and  diagonalized  using  joint  diagonalization, 
the  topic  of  Section  9.3.4.  The  resulting  time-dependent  diagonal  elements  have  been  shown  to 
be  intrinsic  to  the  objects  and  can  be  used,  on  their  own  or  combined  with  other  quantities,  in 
discrimination  processing. 

9.3.4  Joint  diagonalization  for  multi-target  data  pre-processing 

In  real-life  situations  the  targets  of  interest  are  usually  surrounded  by  natural  and  artificial  debris 
with  metallic  content,  including,  for  instance,  the  remains  of  ordnance  that  did  explode.  Thus  it  is 
usually  not  clear  how  many  objects  are  producing  a  given  detected  signal;  all  sensing  methods, 
including  EMI,  are  fraught  with  detection  rates  that  overwhelm  cleanup  efforts  and  hike  their 
cost.  Here  we  introduce  a  data  pre-processing  technique  based  on  joint  diagonalization  (JD)  that 
estimates  the  number  of  targets  present  in  the  field  of  view  of  the  sensor  as  it  takes  a  data  shot, 
and,  in  a  good  number  of  cases,  even  provides  the  capability  to  perform  real-time 
characterization  and  classification  of  the  targets  without  the  need  for  a  forward  model. 

Joint  diagonalization  has  become  an  important  tool  for  signal  processing  and  inverse  problems, 
used  as  part  of  independent  component  analysis  [22],  blind  source  separation  or  BSS  [36], 
common  principal  component  analysis,  and,  more  recently,  kernel-based  nonlinear  BSS  [37].  We 
further  extend  the  applicability  of  the  method  by  using  it  to  detect  and  locate  buried  targets 
without  the  need  for  inversion.  As  we  say  above,  a  variation  of  the  method  can  be  used  to 
extricate  time-dependent  electromagnetic  signatures  from  attitude  information.  Here  we  will 
outline  the  detailed  procedure  as  applied  to  the  TEMTADS  sensor  array,  a  time-domain  device 
with  25  transmitter/receiver  pairs  that  provides  625  measurements  over  Ng  =  123  time  gates  at 
each  sensor  location.  Equation  Section  (Next) 

9.3.5  The  multi- static  response  matrix 

JD  estimates  the  eigenvalues  and  eigenvectors  of  a  square  time-  or  frequency-dependent  multi¬ 
static  response  (MSR)  matrix  synthesized  directly  from  measured  values.  To  construct  the  MSR 
matrices  one  just  has  to  stack  the  625  readings  at  each  time  gate  in  a  25  x  25  array  so  that  each 
column  stands  for  one  of  N,  transmitters  and  each  row  represents  one  of  Nr  receivers: 


S(4): 


Hn 

h12  • 

1 

af 

H2l 

H22  • 

... 

„5S 

,  k  =  l,...,Ng, 

(B.l) 

HNr2  • 

-  HKN,_ 
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where  the  element  HtJ  is  the  field  measured  by  the  /-th  receiver  when  the  j- th  transmitter  is  fired. 
The  second  step  of  the  procedure  is  to  diagonalize  the  123  matrices  at  one  stroke  so  they  all  share 
a  single  set  of  orthonormal  eigenvectors.  In  other  words,  given  the  MSR  matrix  S(4)  at  the  k- th 
time  gate,  we  look  for  a  unitary  matrix  V  such  that  the  products 


Dk=\TS(tk)\  (B.2) 

are  “as  diagonal  as  possible”  (i.e.,  their  off-diagonal  elements  vanish  within  a  preset  tolerance). 
By  diagonalizing  all  the  matrices  simultaneously  we  separate  the  time-dependent  intrinsic 
features  of  the  responding  sources  (and  hence  the  interred  objects),  which  get  encapsulated  in  the 
eigenvalues,  from  the  other  factors — notably  the  location  and  orientation  of  the  target  with 
respect  to  the  sensor — that  influence  the  signal  but  do  not  change  as  the  data  are  being  taken; 
these  get  bundled  into  the  eigenvectors.  (The  fact  that  the  locations  and  orientations  can  be 
dissociated  in  this  way  from  the  electromagnetic  signatures  is  an  upside  of  the  low  frequencies  of 
the  quasistatic  EMI  range,  because  the  relevant  Green  functions  are  time -independent.)  Thus  the 
measured  data  can  be  resolved  as  a  superposition  of  “elemental”  sub-signals,  each  corresponding 
to  an  elementary  dipolar  source,  whose  combination  corresponds  to  the  buried  objects.  Each 
source — and  the  corresponding  field  singularity — can  moreover  be  localized  numerically:  the 
TEMTADS  geometry  is  such  that  the  diagonal  of  the  unprocessed  MSR  matrix  mimics  a  set  of 
monostatic  measurements,  akin  to  those  taken  with  a  handheld  sensor,  which  peak  sharply  when 
there  is  a  target  directly  underneath.  The  maxima  in  the  diagonal  thus  point  to  the 
transmitter/receiver  pairs  closest  to  any  responding  sources.  These  location  estimates  can  be 
grouped  and  correlated  to  the  eigenvalue  distributions  to  estimate  target  locations. 

9.3.6  Interpretation  and  diagonalization  of  the  MSR  matrix 

We  now  proceed  to  express  our  above  considerations  quantitatively.  Initially  we  consider  the 
transmitter  assembly,  which  in  TEMTADS  consists  of  a  set  of  coplanar  square  loops  forming  a 
regular  grid.  The  Biot-Savart  law  gives  the  primary  magnetic  induction  established  at  the 
location  r;  of  the  /-th  source  when  the  j- th  transmitter  antenna  (whose  area  is  <rTx  )  is  excited 

immediately  before  shutoff  by  a  current  If. 


bp: 

ji 


0  j 


■  cr 

\n  Tx 


1  f  dl'x(r ,-r') 


A 


lr(-r'P 


<<Vr 


This  primary  field  induces  in  the  /-th  source  a  dipole  moment  given  by 


(B.3) 


m„  =  U,A,u;  BJ ,  (B.4) 

where  the  Euler  rotation  matrix  U  relates  the  instrument’s  coordinate  axes  to  the  principal  axes 
of  the  source,  and  the  diagonal  polarizability  matrix  A/,  the  only  quantity  intrinsic  to  the  source, 
measures  the  strength  with  which  the  primary  field  induces  a  moment  along  each  of  those  axes. 


According  to  Faraday’s  law,  the  signal  measured  by  a  receiver  coil  is  the  electromotive  force 
given  by  the  negative  of  the  time  derivative  of  the  secondary  magnetic  flux  through  the  coil. 
Since  the  field  at  point  r  of  a  dipole  of  moment  m  placed  at  r0  is  given  by 
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B  =  —V x |  mx 
An 


r-r 


lr-r0l 


^■j,  and  thus  jB-ds  =  -m--^-(|)dl 


I r  _  ro  I 


by  straightforward  application  of  Stokes’s  theorem,  one  obtains  that  the  signal  sampled  at  time  4 
by  the  z-th  receiver  (of  area  crRx  )  when  the  /- th  source  is  excited  by  the  j- th  transmitter  is 


H  ,j  (4  )CrRxi  ^Tx  ,  I  j 


M0 


1  r  d\'i  x  (r'  -  r( ) 


■  <j) 


4^  KX'A7RXj/Xj  I  r'  -  r,  I3 
=  g;>RXi-[UA((4)Ur]-g^T  /., 


hi ,,  (4 )  =  g(ccrRx  •  hi (4 ) 


(B.5) 


where  a  dot  over  a  variable  indicates  its  time  derivative.  In  equations  (B.3)  and  (B.5)  the  line 
element  d\'  lies  on  the  x-y  plane,  and  as  a  consequence  the  Green  functions  are  similar  in 
structure  to  those  of  the  simple  model  presented  in  Section  2.2.  Note  that  we  have  included  the 
exciting  current  Ij  and  the  transmitter  and  receiver  areas  in  the  definition  of  the  signal;  we  have 
explicit  knowledge  of  these  quantities  and  can  factor  them  out.  If  only  the  l- th  source  is 
illuminated,  we  construct  the  MSR  matrix  for  the  complete  transmitter/receiver  array  by  tiling  Nr 
x  Nt  instances  of  the  expression  (B.5): 


S  =  GscU;A,U[(Gpr)r,  (B.6) 

where  the  primary  (or  transmitter)  dyad  Gpr  is  of  size  Ntx  3,  the  secondary  (or  receiver)  dyad  Gsc 
is  of  size  Nr  x  3,  and  the  response  matrix  UA/Ur  is  3  x  3.  When  there  is  more  than  one  source 
present,  the  MSR  matrix  of  equation  (B.6)  is  readily  generalized: 


"uAuf 
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(B.7) 


where  we  see  that  the  features  intrinsic  to  the  targets  can  be  separated  formally  from  the 
particulars  of  the  measurement — that  is,  from  the  geometry  and  dimensions  of  the  sensor  and  the 
sensor-target  attitude.  The  array  S  has  size  Nr  x  N,  and  is  square  if  Nr  =  Nt,  as  is  the  case  with 
TEMTADS.  This  allows  us  to  diagonalize  the  matrix  but  does  not  suffice  to  guarantee  that  the 
extracted  information  is  useful — i.e.,  that  the  eigenvalues  and  eigenvectors  are  real,  and  that  the 
latter  are  orthonormal.  For  that  to  hold  we  must  have  a  real,  symmetric  matrix,  which  requires 

G*c  =  Gpr  =  G; .  This  cannot  be  rigorously  true,  because  the  receivers  cannot  coincide  exactly 

with  the  transmitters,  but  holds  approximately  for  TEMTADS  if  we  factor  the  exciting  current 
and  the  coil  areas  out  of  S,  as  we  did  in  equation  (B.5).  The  diagonalization  we  perform  is  thus  a 
particular  case  of  a  singular  value  decomposition  (SVD),  and  in  what  follows  we  use 
“diagonalization”  as  shorthand  for  “SVD  of  a  symmetric  matrix.” 
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The  decomposition  (B.7)  exhibits  the  actual  polarizability  elements  but  is  not  directly  available 
to  us  because  the  Green  tensors  are  not  orthogonal.  To  see  what  we  do  get  when  we  diagonalize 
S  we  can  perform  the  SVD  on  G: 

S  =  GUAUrGr  =  w[zVrUAUrVz]  Wr  =  WZA7/W7  =  YAYr  (B.8) 

In  the  intermediate  step  we  have  used  the  fact  that  the  matrix  within  the  brackets  is  real  and 
symmetric  and  thus  has  a  purely  real  eigendecomposition.  Result  (B.8)  shows  that  the  eigenvalue 
matrix  A,  though  time-dependent,  is  not  solely  composed  of  source  responses,  but  also  contains 
location  and  orientation  information  extracted  from  the  Green  tensors.  The  eigenvectors, 
likewise,  include  information  from  both  the  polarizabilities  and  the  measurement  particulars. 

We  also  see  in  the  decomposition  (B.8)  that  S  contains  an  unknown  “hidden  dimension” — 3 N, 
where  N  is  the  number  of  sources — in  the  size  of  the  block-diagonal  response  matrix.  Numerical 
diagonalization  (or,  in  general,  the  SVD)  of  S  will  impose  this  middle  dimension  to  be  Nr  =  Nt. 
Ideally,  the  method  should  be  able  to  resolve  up  to  [ Nr  /  3j  responding  sources,  or  eight  for 

TEMTADS,  but  the  actual  number  is  lower.  For  one,  the  procedure  will  resolve  targets  only 
when  they  are  spatially  separated:  two  distinct  dipoles  sharing  one  location  decrease  the  rank  of 
the  G  matrices,  and  hence  of  S,  by  3.  In  any  case,  diagonalization  of  S  can  again  let  us  estimate 
the  number  of  targets  illuminated  by  the  sensor;  since  the  only  time-dependent  quantities  are  the 
intrinsic  polarizabilities  of  the  sources,  we  expect  the  additional  information  provided  by  the 
time  decay  of  the  eigenvalues  to  be  useful  for  classification. 

The  development  outlined  above  corresponds  to  each  time  gate  taken  separately.  To  make  sense 
of  the  time-dependent  information  we  have  to  find  a  way  to  “follow”  each  of  the  eigenvalues  as 
the  signal  decays.  (A  similar  process  must  be  carried  out  when  using  the  dipole  model  for 
inversion.)  One  could  in  principle  diagonalize  the  MSR  matrix  at  each  time  channel,  and  the 
eigenvectors,  which  depend  only  on  geometry  and  pose,  should  stay  constant;  however,  it  is  not 
possible  to  know  a  priori  the  order  in  which  the  eigenvalues  will  be  given  by  the  diagonalization; 
this  fact — not  to  mention  noise  and  experimental  uncertainty — makes  it  inevitable  to  have  to 
disentangle  the  tensor  elements  by  hand,  which  is  easily  done  wrong.  Instead,  we  explicitly  look 
for  an  orthogonal  matrix  of  eigenvectors  that  diagonalizes  all  the  MSR  matrices  simultaneously. 
The  procedure  we  employ  is  a  generalization  of  the  method  for  single  matrices,  and  is  well- 
known;  it  is  sketched  in  next  Section. 

9. 3. 7  Algorithm  for  joint  diagonalization 

The  joint  diagonalization  algorithm  we  use  [36,  38-39]  is  a  generalization  of  Jacobi’s  procedure 
to  find  the  eigenvalues  of  a  single  matrix.  Formally  we  set  out  to  solve  the  optimization  problem 


min 

V 

s.t. 


“  q= 1  ‘*j 


=  /, 


(B.9) 


which  we  accomplish  by  making  repeated  Givens -Jacobi  similarity  transformations  designed  to 
gradually  accumulate  the  “content”  of  the  matrices  on  their  diagonals  until  a  certain  tolerance 
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level  is  reached.  The  transformations  are  of  the  form  A(t  )  —>  A'(t  )  =  V  A(t  )VT ,  with  the  matrix 

v  q  7  q  rs  q  rs  7 

Vn  being  the  identity  but  with  the  four  elements  Vrr,  Vrs,  Vsr,  and  Vss  replaced  by  the  two- 
dimensional  rotation  array 
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+  n 


(B.10) 


where 


n  (t  )-  a  (t  -[a  ( t  )  + a  ( t  )]2 1,  (B.ll) 

rs  7  J  f -  rr  x  q y  ss  v  q ' J  L  rs v  qy  sr  v  q 7  A  y  v  y 

q 

f  =2 Y[a  ( t)-a  (. t  )][a  ( t)  +  a  (t  )] .  (B.12) 

J  rs  L  rr  v  q 7  ss  v  q rs  v  q '  sr  x  q/J  v  7 

q 

The  indices  are  swept  systematically,  and  the  procedure  is  repeated  until  convergence  is  reached. 
The  computational  burden  is  equivalent  to  that  of  diagonalizing  the  matrices  one  by  one.  The 
resulting  eigenvalues  and  eigenvectors  are  all  real  because  all  the  MSR  matrices  are  symmetric. 


43 


April  2012 


